﻿//pore-scale study on reactive mixing of miscible with viscous fingering in porous media�����и�д
//MRT-LBM for N-S equation ; LKS for ADR equation.
//simulation of the


#define _CRT_SECURE_NO_WARNINGS
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define Nx 1024         //numbers of cells in the x-direction
#define Ny 127          //numbers of cells in the y-direction
#define Nx1 (Nx+1)      //numbers of points in the x-direction
#define Ny1 (Ny+1)      //numbers of points in the y-direction
#define Q 9             //numbers of discrete velocities
#define rho0 1.0        
#define H 32.0            
#define L 256.0           
#define Pe_a 200       
#define Pe_b 200      
#define Pe_c 200        
#define Pe_d 2000        
#define niu_a 0.005
#define niu_b 0.5
#define niu_c 0.05
//#define niu_d 0.0005
#define tau 1.0
#define U0 0.005        
#define step 80000      
#define step_out 20000	//output step


double rho[Nx1][Ny1], u[Nx1][Ny1][2], u0[Nx1][Ny1][2], f[Nx1][Ny1][Q], f_post[Nx1][Ny1][Q];
double g_a[Nx1][Ny1][Q],g_ap[Nx1][Ny1][Q],g_b[Nx1][Ny1][Q],g_bp[Nx1][Ny1][Q],g_c[Nx1][Ny1][Q],g_cp[Nx1][Ny1][Q],g_d[Nx1][Ny1][Q],g_dp[Nx1][Ny1][Q];
double r[Nx1][Ny1][Q],R[Nx1][Ny1][Q],rs[Nx1][Ny1][Q],RS[Nx1][Ny1][Q];
double c_a[Nx1][Ny1],c_b[Nx1][Ny1],c_c[Nx1][Ny1],c_d[Nx1][Ny1],tau1[Nx1][Ny1],niu[Nx1][Ny1],viu[Nx1][Ny1];
double dx,dy,dt,maxx,cc,cs,Dm_a,Dm_b,Dm_c,Dm_d,A_a,A_b,A_c,A_d;
double s[Q];
int L_d;                
char data[20],matlab[20],data_a[20],data_b[20],data_c[20],data_d[20];
double F[2]={0,0};
double D[Q] = { 9,36,36,6,12,6,12,4,4 };
double w[Q] = { 4.0 / 9,1.0 / 9,1.0 / 9,1.0 / 9,1.0 / 9,1.0 / 36,1.0 / 36,1.0 / 36,1.0 / 36 };
int e[Q][2] = { {0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1} };


//===============================================================================================================================
double feq(double RHO, double u[2], int k)              
{
	double cu, U2;
	cu = e[k][0] * u[0] + e[k][1] * u[1];
	U2 = u[0] * u[0] + u[1] * u[1];
	return w[k] * (RHO + 3.0*cu + 4.5*cu*cu - 1.5*U2);
}

double reaction(int i, int j, int k)                    
{
	double uu, a, m;
	a = 0.05;
	uu = a * c_a[i][j] * c_b[i][j];
	m = w[k] * uu*(1 + (e[k][0] * u[i][j][0] + e[k][1] * u[i][j][1]) / cs / cs / tau * (tau - 0.5));
	return m;
}

double reactions(int i, int j, int k)
{
	double uu, k_n,k_p, m,c_n,c_s;                      
	c_n = 0.15;
	c_s = 0.10;
	if (c_c[i][j] > c_n)
		k_n = 0.0;
	else
		k_n = 0.0;
	if (c_c[i][j] > c_s)
		k_p = 0.0;
	else
		k_p = 0.0;
	uu = k_n * c_c[i][j] * (c_c[i][j] - c_n) + k_p * c_c[i][j] * c_d[i][j] * (c_c[i][j] - c_s);
	m = w[k] * uu *(1 + (e[k][0] * u[i][j][0] + e[k][1] * u[i][j][1]) / cs / cs * (tau - 0.5) / tau);
	return m;
}

double geq_a(int i,int j,int k)                           //����A�ķֲ�������ʼ��
{
	double cu, U2;
	cu = e[k][0] * u[i][j][0] + e[k][1] * u[i][j][1];
	U2 = u[i][j][0] * u[i][j][0] + u[i][j][1] * u[i][j][1];
	return w[k] * c_a[i][j] * (1.0 + 3.0* cu + 4.5*cu*cu - 1.5*U2);
}

double geq_b(int i,int j,int k)                          
{
	double cu, U2;
	cu = e[k][0] * u[i][j][0] + e[k][1] * u[i][j][1];
	U2 = u[i][j][0] * u[i][j][0] + u[i][j][1] * u[i][j][1];
	return w[k] * c_b[i][j] * (1.0 + 3.0* cu + 4.5*cu*cu - 1.5*U2);
}

double geq_c(int i,int j,int k)                          
{
	double cu, U2;
	cu = e[k][0] * u[i][j][0] + e[k][1] * u[i][j][1];
	U2 = u[i][j][0] * u[i][j][0] + u[i][j][1] * u[i][j][1];
	return w[k] * c_c[i][j] * (1.0 + 3.0* cu + 4.5*cu*cu - 1.5*U2);
}

double geq_d(int i,int j,int k)                          
{
	double cu, U2;
	cu = e[k][0] * u[i][j][0] + e[k][1] * u[i][j][1];
	U2 = u[i][j][0] * u[i][j][0] + u[i][j][1] * u[i][j][1];
	return w[k] * c_d[i][j] * (1.0 + 3.0* cu + 4.5*cu*cu - 1.5*U2);
}

double geq_al(int i, int j, int k)                       
{
	int n;
	double m1, m2,x;
	double cu,uu;
	cu = e[k][0] * u[i][j][0] + e[k][1] * u[i][j][1];
	uu = u[i][j][0] * u[i][j][0] + u[i][j][1] * u[i][j][1];
	m1 = 0.0; m2 = 0.0;
	for (n = 0; n < Q; n++)
	{
		m1 = m1 + e[n][0] * g_a[i][j][n] + tau* dt * e[n][0] * R[i][j][n];
		m2 = m2 + e[n][1] * g_a[i][j][n] + tau* dt * e[n][1] * R[i][j][n];
	}
	x = w[k] * c_a[i][j] * (1.0 + cu / cs / cs) + w[k] * A_a / (cs*cs) / (A_a - tau)*(e[k][0] * (m1 - c_a[i][j] * u[i][j][0])
		+ e[k][1] * (m2 - c_a[i][j] * u[i][j][1]));
	return x;
}

double geq_bl(int i, int j, int k)                        
{
	int n;
	double m1, m2,x;
	double cu,uu;
	cu = e[k][0] * u[i][j][0] + e[k][1] * u[i][j][1];
	uu = u[i][j][0] * u[i][j][0] + u[i][j][1] * u[i][j][1];
	m1 = 0.0; m2 = 0.0;
	for (n = 0; n < Q; n++)
	{
		m1 = m1 + e[n][0] * g_b[i][j][n] + tau* dt * e[n][0] * R[i][j][n];
		m2 = m2 + e[n][1] * g_b[i][j][n] + tau* dt * e[n][1] * R[i][j][n];
	}
	x = w[k] * c_b[i][j] * (1.0 + cu / cs / cs) + w[k] * A_b / (cs*cs) / (A_b - tau)*(e[k][0] * (m1 - c_b[i][j] * u[i][j][0])
		+ e[k][1] * (m2 - c_b[i][j] * u[i][j][1]));
	return x;
}

double geq_cl(int i, int j, int k)                        
{
	int n;
	double m1, m2,x;
	double cu,uu;
	cu = e[k][0] * u[i][j][0] + e[k][1] * u[i][j][1];
	uu = u[i][j][0] * u[i][j][0] + u[i][j][1] * u[i][j][1];
	m1 = 0.0; m2 = 0.0;
	for (n = 0; n < Q; n++)
	{
		m1 = m1 + e[n][0] * g_c[i][j][n] - tau* dt * e[n][0] * r[i][j][n]+ tau*dt*e[n][0]*rs[i][j][n];
		m2 = m2 + e[n][1] * g_c[i][j][n] - tau* dt * e[n][1] * r[i][j][n]+ tau*dt*e[n][1]*rs[i][j][n];
	}
	x = w[k] * c_c[i][j] * (1.0 + cu / cs / cs) + w[k] * A_c / (cs*cs) / (A_c - tau)*(e[k][0] * (m1 - c_c[i][j] * u[i][j][0])
		+ e[k][1] * (m2 - c_c[i][j] * u[i][j][1]));
	return x;
}

double geq_dl(int i, int j, int k)                       
{
	int n;
	double m1, m2,x;
	double cu,uu;
	cu = e[k][0] * u[i][j][0] + e[k][1] * u[i][j][1];
	uu = u[i][j][0] * u[i][j][0] + u[i][j][1] * u[i][j][1];
	m1 = 0.0; m2 = 0.0;
	for (n = 0; n < Q; n++)
	{
		m1 = m1 + e[n][0] * g_d[i][j][n] - tau* dt * e[n][0] * rs[i][j][n];
		m2 = m2 + e[n][1] * g_d[i][j][n] - tau* dt * e[n][1] * rs[i][j][n];
	}
	x = w[k] * c_d[i][j] * (1.0 + cu / cs / cs) + w[k] * A_d / (cs*cs) / (A_d - tau)*(e[k][0] * (m1 - c_d[i][j] * u[i][j][0])
		+ e[k][1] * (m2 - c_d[i][j] * u[i][j][1]));
	return x;
}


//===============================================================================================================================
void Init_Eq()                  
{
	int i, j, k;
	for (i = 0; i <= Nx; i++) for (j = 0; j <= Ny; j++)
	{
		c_a[i][j] = 0.0;
		c_b[i][j] = 1.0;
		c_c[i][j] = 0.0;
		c_d[i][j] = 0.0;
		rho[i][j] = rho0;
		u[i][j][0] = 0.0;
		u[i][j][1] = 0.0;

		for (k = 0; k < Q; k++)
		{
			f[i][j][k] = feq(rho[i][j], u[i][j], k);
			g_a[i][j][k] = geq_a(i,j , k);
			g_b[i][j][k] = geq_b(i,j , k);
			g_c[i][j][k] = geq_c(i,j , k);
			g_d[i][j][k] = geq_d(i,j , k);
			R[i][j][k] = reaction(i,j,k);
			RS[i][j][k]=reactions(i,j,k);
		}
	}
	for(i=0;i<=0.1*Nx;i++)
		for (j = 0; j < Ny1; j++)
		{
			c_a[i][j] = 1.0;
			c_b[i][j] = 0.0;
			c_c[i][j] = 0.0;
			c_d[i][j] = 0.0;
			u[i][j][0] = U0;
			u[i][j][1] = 0.0;
			rho[i][j] = rho0;
			for (k = 0; k < Q; k++)
			{
				f[i][j][k] = feq(rho[i][j], u[i][j], k);
                g_a[i][j][k] = geq_a(i,j , k);
                g_b[i][j][k] = geq_b(i,j , k);
                g_c[i][j][k] = geq_c(i,j , k);
                g_d[i][j][k] = geq_d(i,j , k);
                R[i][j][k] = reaction(i,j,k);
                RS[i][j][k]= reactions(i,j,k);
			}
		}
}


//===============================================================================================================================
void coll_BGK()
{
	double GEQ_a,GEQ_b,GEQ_c,GEQ_d;
	int i, j, k;
	for(i=0;i<Nx1;i++)
		for(j=0;j<Ny1;j++)
			for (k = 0; k < Q; k++)
			{
				r[i][j][k]= reaction(i, j, k);
				rs[i][j][k] = reactions(i,j,k);

				GEQ_a = geq_al(i, j, k);
				GEQ_b = geq_bl(i, j, k);
				GEQ_c = geq_cl(i, j, k);
				GEQ_d = geq_dl(i, j, k);

				g_ap[i][j][k] = g_a[i][j][k]-(g_a[i][j][k]-GEQ_a)/tau - dt * r[i][j][k] - dt / 2 * (r[i][j][k] - R[i][j][k]);
				g_bp[i][j][k] = g_b[i][j][k]-(g_b[i][j][k]-GEQ_b)/tau - dt * r[i][j][k] - dt / 2 * (r[i][j][k] - R[i][j][k]);
				g_cp[i][j][k] = g_c[i][j][k]-(g_c[i][j][k]-GEQ_c)/tau + dt * (r[i][j][k]-rs[i][j][k]) + dt / 2 * (r[i][j][k] - R[i][j][k]) -
                    dt / 2 * (rs[i][j][k]-RS[i][j][k]);
				g_dp[i][j][k] = g_d[i][j][k]-(g_d[i][j][k]-GEQ_d)/tau + dt * rs[i][j][k] + dt / 2 * (rs[i][j][k] - RS[i][j][k]);

				R[i][j][k]=r[i][j][k];
				RS[i][j][k]=rs[i][j][k];             
			}
}


//===============================================================================================================================
double meq(double P, double u[2], int k)
{
	double x;
	switch (k)
	{
	case 0: {x =P; break; }
	case 1: {x =(-2*P+3*P*(u[1]*u[1]+u[0]*u[0])); break; }
	case 2: {x =(P-3*P*(u[0]*u[0]+u[1]*u[1])); break; }
	case 3: {x =P* u[0]; break; }
	case 4: {x = -P* u[0]; break; }
	case 5: {x =P* u[1]; break; }
	case 6: {x = -P* u[1]; break; }
	case 7: {x = P* (u[0]*u[0] - u[1]* u[1]); break; }
	case 8: {x = P* u[0]*u[1]; break; }
	default: x = 0;
	}
	return x;
}

double force(double u[2], int k)
{
	double x;
	switch (k)
	{
		case 0: {x=0; break;}
		case 1: {x=6*(u[0]*F[0]+u[1]*F[1]); break; }
		case 2: {x=-6*(u[0]*F[0]+u[1]*F[1]); break; }
		case 3: {x=F[0]; break; }
		case 4: {x=-F[0]; break; }
		case 5: {x=F[1]; break; }
		case 6: {x=-F[1]; break; }
		case 7: {x=2*(u[0]*F[0]-u[1]*F[1]); break; }
		case 8: {x=u[0]*F[1]+u[1]*F[0]; break; }
		default: x=0;
	}
	return x;
}

void Coll_MRT()
{
	int j, i, k;
	double MEQ,forcing;
	double m[Q];
	for (i= 0; i<= Nx; i++) for (j = 0; j <= Ny; j++)
	{
		// Transformation from velocity space to moment space:
		m[0] = f[i][j][0] + f[i][j][1] + f[i][j][2] + f[i][j][3] + f[i][j][4] + f[i][j][5] + f[i][j][6] + f[i][j][7] + f[i][j][8];
		m[1] = -4 * f[i][j][0] - f[i][j][1] - f[i][j][2] - f[i][j][3] - f[i][j][4] + 2 * (f[i][j][5] + f[i][j][6] + f[i][j][7] + f[i][j][8]);
        m[2] = 4 * f[i][j][0] - 2 * (f[i][j][1] + f[i][j][2] + f[i][j][3] + f[i][j][4]) + f[i][j][5] + f[i][j][6] + f[i][j][7] + f[i][j][8];
        m[3] = f[i][j][1] - f[i][j][3] + f[i][j][5] - f[i][j][6] -f[i][j][7] + f[i][j][8];
        m[4] = -2 * (f[i][j][1] - f[i][j][3]) + f[i][j][5] - f[i][j][6] - f[i][j][7] + f[i][j][8];
        m[5] = f[i][j][2] - f[i][j][4] + f[i][j][5] + f[i][j][6] - f[i][j][7] - f[i][j][8];
        m[6] = -2 * (f[i][j][2] - f[i][j][4]) + f[i][j][5] + f[i][j][6] -f[i][j][7] - f[i][j][8];
        m[7] = f[i][j][1] - f[i][j][2] + f[i][j][3] - f[i][j][4];
		m[8] = f[i][j][5] - f[i][j][6] + f[i][j][7] - f[i][j][8];

		niu[i][j] = niu_a * exp(c_b[i][j] * log(niu_b/niu_a)+c_c[i][j]*log(niu_c/niu_a));
		viu[i][j] = niu[i][j] / rho[i][j];
		tau1[i][j] = 3 * viu[i][j]/dt + 0.5;

		s[7] = s[8] = 1.0 / tau1[i][j]; s[0] = s[3] = s[5] = 0.0; s[4] = s[6] = 8.0 * (2.0 - s[7]) / (8.0 - s[7]);
		s[1] = 1.0/tau1[i][j]; s[2] = 1.0/tau1[i][j];

		for (k = 0; k < Q; k++)
		{
			MEQ = meq(rho[i][j], u[i][j], k);
			forcing = force(u[i][j],k);
			m[k] = m[k] - s[k] * (m[k] - MEQ) + dt * forcing;           // relaxation
			m[k] /= D[k];                                // rescaling
		}

		f_post[i][j][0] = m[0] - 4 * (m[1] - m[2]);
		f_post[i][j][1] = m[0] - m[1] - 2 * (m[2] + m[4]) + m[3] + m[7];
		f_post[i][j][2] = m[0] - m[1] - 2 * (m[2] + m[6]) + m[5] - m[7];
		f_post[i][j][3] = m[0] - m[1] - 2 * (m[2] - m[4]) - m[3] + m[7];
		f_post[i][j][4] = m[0] - m[1] - 2 * (m[2] - m[6]) - m[5] - m[7];
		f_post[i][j][5] = m[0] + m[1] + m[1] + m[2] + m[3] + m[4] + m[5] + m[6]+ m[8];
		f_post[i][j][6] = m[0] + m[1] + m[1] + m[2] - m[3] - m[4] + m[5] + m[6]- m[8];
		f_post[i][j][7] = m[0] + m[1] + m[1] + m[2] - m[3] - m[4] - m[5] - m[6]+ m[8];
		f_post[i][j][8] = m[0] + m[1] + m[1] + m[2] + m[3] + m[4] - m[5] - m[6]- m[8];
	}
}


//===============================================================================================================================
void Streaming()
{
	int i, j, jd, id, k;
	for (i = 0; i <= Nx; i++) for (j = 0; j <= Ny; j++) for (k = 0; k < Q; k++)
	{
		id = i - e[k][0]; jd = j - e[k][1];                 // upwind node
		if (jd >= 0 && jd <= Ny && id >= 0 && id <= Nx)     // fluid node
		{
			f[i][j][k] = f_post[id][jd][k];                 // streaming
			g_a[i][j][k] = g_ap[id][jd][k];
			g_b[i][j][k] = g_bp[id][jd][k];
			g_c[i][j][k] = g_cp[id][jd][k];
			g_d[i][j][k] = g_dp[id][jd][k];
		}
	}
}


//===============================================================================================================================
void Boundary()                                             
{
	int i, j,k;

	for (i = 0; i <=Nx; i++)                                
	{
		f[i][0][2] = f_post[i][0][4];
		f[i][0][5] = f_post[i][0][7];
		f[i][0][6] = f_post[i][0][8];

		g_a[i][0][2] = g_ap[i][0][4];
		g_a[i][0][5] = g_ap[i][0][7];
		g_a[i][0][6] = g_ap[i][0][8];

		g_b[i][0][2] = g_bp[i][0][4];
		g_b[i][0][5] = g_bp[i][0][7];
		g_b[i][0][6] = g_bp[i][0][8];

		g_c[i][0][2] = g_cp[i][0][4];
		g_c[i][0][5] = g_cp[i][0][7];
		g_c[i][0][6] = g_cp[i][0][8];

		g_d[i][0][2] = g_dp[i][0][4];
		g_d[i][0][5] = g_dp[i][0][7];
		g_d[i][0][6] = g_dp[i][0][8];
	}

	for (i = 0; i <= Nx; i++)                               
	{
        f[i][Ny][4] = f_post[i][Ny][2];
		f[i][Ny][7] = f_post[i][Ny][5];
		f[i][Ny][8] = f_post[i][Ny][6];

		g_a[i][Ny][4] = g_ap[i][Ny][2];
		g_a[i][Ny][7] = g_ap[i][Ny][5];
		g_a[i][Ny][8] = g_ap[i][Ny][6];

		g_b[i][Ny][4] = g_bp[i][Ny][2];
		g_b[i][Ny][7] = g_bp[i][Ny][5];
		g_b[i][Ny][8] = g_bp[i][Ny][6];

		g_c[i][Ny][4] = g_cp[i][Ny][2];
		g_c[i][Ny][7] = g_cp[i][Ny][5];
		g_c[i][Ny][8] = g_cp[i][Ny][6];

		g_d[i][Ny][4] = g_dp[i][Ny][2];
		g_d[i][Ny][7] = g_dp[i][Ny][5];
		g_d[i][Ny][8] = g_dp[i][Ny][6];
	}

	for (i = 1; i < Nx; i++) for (j = 0; j <= Ny; j++)      
	{
		rho[i][j] = f[i][j][0] + f[i][j][1] + f[i][j][2] + f[i][j][3] + f[i][j][4] + f[i][j][5] + f[i][j][6] + f[i][j][7] + f[i][j][8];
		u[i][j][0] = (f[i][j][1] + f[i][j][5] + f[i][j][8] - f[i][j][3] - f[i][j][6] - f[i][j][7]) / rho[i][j];
		u[i][j][1] = (f[i][j][2] + f[i][j][5] + f[i][j][6] - f[i][j][4] - f[i][j][7] - f[i][j][8]) / rho[i][j];
		c_a[i][j] = g_a[i][j][0] + g_a[i][j][1] + g_a[i][j][2] + g_a[i][j][3] + g_a[i][j][4] + g_a[i][j][5] + g_a[i][j][6] + g_a[i][j][7] + g_a[i][j][8];
		c_b[i][j] = g_b[i][j][0] + g_b[i][j][1] + g_b[i][j][2] + g_b[i][j][3] + g_b[i][j][4] + g_b[i][j][5] + g_b[i][j][6] + g_b[i][j][7] + g_b[i][j][8];
		c_c[i][j] = g_c[i][j][0] + g_c[i][j][1] + g_c[i][j][2] + g_c[i][j][3] + g_c[i][j][4] + g_c[i][j][5] + g_c[i][j][6] + g_c[i][j][7] + g_c[i][j][8];
		c_d[i][j] = g_d[i][j][0] + g_d[i][j][1] + g_d[i][j][2] + g_d[i][j][3] + g_d[i][j][4] + g_d[i][j][5] + g_d[i][j][6] + g_d[i][j][7] + g_d[i][j][8];
	}

	for (j = 0; j <= Ny; j++)                               
	{
		u[0][j][0] = U0;
		u[0][j][1] = 0.0;
		c_a[0][j] = 1.0;
		c_b[0][j] = 0.0;
		c_c[0][j] = 0.0;
		c_d[0][j] = 0.0;

		for (k = 0; k < Q; k++)
		{

			f[0][j][1] = feq(rho[1][j], u[0][j], 1) + f[1][j][1] - feq(rho[1][j], u[1][j], 1);
			f[0][j][5] = feq(rho[1][j], u[0][j], 5) + f[1][j][5] - feq(rho[1][j], u[1][j], 5);
			f[0][j][8] = feq(rho[1][j], u[0][j], 8) + f[1][j][8] - feq(rho[1][j], u[1][j], 8);


			g_a[0][j][1] = geq_al(0, j, 1) + g_a[1][j][1] - geq_al(1, j, 1);
			g_a[0][j][5] = geq_al(0, j, 5) + g_a[1][j][5] - geq_al(1, j, 5);
			g_a[0][j][8] = geq_al(0, j, 8) + g_a[1][j][8] - geq_al(1, j, 8);

			g_b[0][j][1] = geq_bl(0, j, 1) + g_b[1][j][1] - geq_bl(1, j, 1);
			g_b[0][j][5] = geq_bl(0, j, 5) + g_b[1][j][5] - geq_bl(1, j, 5);
			g_b[0][j][8] = geq_bl(0, j, 8) + g_b[1][j][8] - geq_bl(1, j, 8);

			g_c[0][j][1] = geq_cl(0, j, 1) + g_c[1][j][1] - geq_cl(1, j, 1);
			g_c[0][j][5] = geq_cl(0, j, 5) + g_c[1][j][5] - geq_cl(1, j, 5);
			g_c[0][j][8] = geq_cl(0, j, 8) + g_c[1][j][8] - geq_cl(1, j, 8);

			g_d[0][j][1] = geq_dl(0, j, 1) + g_d[1][j][1] - geq_dl(1, j, 1);
			g_d[0][j][5] = geq_dl(0, j, 5) + g_d[1][j][5] - geq_dl(1, j, 5);
			g_d[0][j][8] = geq_dl(0, j, 8) + g_d[1][j][8] - geq_dl(1, j, 8);
		}

		rho[0][j] = f[0][j][0] + f[0][j][1] + f[0][j][2] + f[0][j][3] + f[0][j][4] + f[0][j][5] + f[0][j][6] + f[0][j][7] + f[0][j][8];
	}

	for (j = 0; j <=Ny; j++)                                
	{
		maxx = 0.0;
		if (u[Nx - 1][j][0] > maxx)
		{
			maxx = u[Nx-1][j][0];
		}
	}

	for (j = 0; j < Ny1; j++)
	{
		for (k = 0; k < Q; k++)
		{
			f[Nx][j][k] = (f[Nx][j][k] + maxx * f[Nx - 1][j][k]) / (1 + maxx);

			g_a[Nx][j][k] = (g_a[Nx][j][k] + maxx * g_a[Nx - 1][j][k]) / (1 + maxx);
			g_b[Nx][j][k] = (g_b[Nx][j][k] + maxx * g_b[Nx - 1][j][k]) / (1 + maxx);
			g_c[Nx][j][k] = (g_c[Nx][j][k] + maxx * g_c[Nx - 1][j][k]) / (1 + maxx);
			g_d[Nx][j][k] = (g_d[Nx][j][k] + maxx * g_d[Nx - 1][j][k]) / (1 + maxx);
		}
	}

	i = Nx;
	for (j = 0; j <= Ny; j++)
	{
			rho[i][j] = f[i][j][0] + f[i][j][1] + f[i][j][2] + f[i][j][3] + f[i][j][4] + f[i][j][5] + f[i][j][6] + f[i][j][7] + f[i][j][8];
			u[i][j][0] = (f[i][j][1] + f[i][j][5] + f[i][j][8] - f[i][j][3] - f[i][j][6] - f[i][j][7]) / rho[i][j];
			u[i][j][1] = (f[i][j][2] + f[i][j][5] + f[i][j][6] - f[i][j][4] - f[i][j][7] - f[i][j][8]) / rho[i][j];
			c_a[i][j] = g_a[i][j][0] + g_a[i][j][1] + g_a[i][j][2] + g_a[i][j][3] + g_a[i][j][4] + g_a[i][j][5] + g_a[i][j][6] + g_a[i][j][7] + g_a[i][j][8];
            c_b[i][j] = g_b[i][j][0] + g_b[i][j][1] + g_b[i][j][2] + g_b[i][j][3] + g_b[i][j][4] + g_b[i][j][5] + g_b[i][j][6] + g_b[i][j][7] + g_b[i][j][8];
            c_c[i][j] = g_c[i][j][0] + g_c[i][j][1] + g_c[i][j][2] + g_c[i][j][3] + g_c[i][j][4] + g_c[i][j][5] + g_c[i][j][6] + g_c[i][j][7] + g_c[i][j][8];
            c_d[i][j] = g_d[i][j][0] + g_d[i][j][1] + g_d[i][j][2] + g_d[i][j][3] + g_d[i][j][4] + g_d[i][j][5] + g_d[i][j][6] + g_d[i][j][7] + g_d[i][j][8];
	}

}


//================================================================================================================================
void data_output()
{
	int i, j;
	double sum;
	FILE *fp,*fp1,*fp2,*fp3,*fp4,*fp5;

	fp = fopen(matlab, "w+");               //use for matlab
	for (j = 0; j <= Ny; j++) for (i = 0; i <= Nx; i++)
		fprintf(fp, "%d %d %e %e %e %e %e %e \n", i, j, u[i][j][0], u[i][j][1], c_a[i][j],c_b[i][j],c_c[i][j],c_d[i][j]);
	fclose(fp);

	fp1 = fopen(data, "w+");
	fprintf(fp1, "title='Poiseluille Flow'\n");
	fprintf(fp1, "VARIABLES=\"X\",\"Y\",\"U\",\"V\",\"Ca\",\"Cb\",\"Cd\",\"Cds\"\n");
	fprintf(fp1, "ZONE T=\"BOX\",I=%d,J=%d,F=POINT\n", Nx1, Ny1);
	for (j = 0; j <= Ny; j++) for (i = 0; i <= Nx; i++)
		fprintf(fp1, "%d %d %e %e %e %e %e %e \n", i, j, u[i][j][0], u[i][j][1], c_a[i][j], c_b[i][j], c_c[i][j], c_d[i][j]);
	fclose(fp1);

	fp2 = fopen(data_b, "w+");
	for (i = 0; i <= Nx; i++)
	{
		sum = 0;
		for (j = 0; j <= Ny; j++)
			sum += c_b[i][j];

		fprintf(fp2, "%d %e \n", i, sum / Ny1);

	}
	fclose(fp2);

	fp3 = fopen(data_a, "w+");
	L_d = 0;
	for (i = 0; i <= Nx; i++)
	{
		sum = 0;
		for (j = 0; j <= Ny; j++)
			sum += c_a[i][j];
		sum = sum / Ny1;
		if (sum >= 0.01&&sum <= 0.99)
			L_d = L_d + 1;
		fprintf(fp3, "%d %e \n", i, sum );

	}
	fclose(fp3);

	fp4 = fopen(data_c, "w+");
	for (i = 0; i <= Nx; i++)
	{
		sum = 0;
		for (j = 0; j <= Ny; j++)
			sum += c_c[i][j];

		fprintf(fp4, "%d %e \n", i, sum / Ny1);

	}
	fclose(fp4);

	fp5 = fopen(data_d, "w+");
	for (i = 0; i <= Nx; i++)
	{
		sum = 0;
		for (j = 0; j <= Ny; j++)
			sum += c_d[i][j];

		fprintf(fp5, "%d %e \n", i, sum / Ny1);

	}
	fclose(fp5);
}


//================================================================================================================================
int main()
{
	int k;
	double err;

	dx = H / (Ny+1);
	dt = dx;
	cc = dx / dt;
	cs = cc / sqrt(3.0);

	Dm_a = U0 * H / Pe_a;
	Dm_b = U0 * H / Pe_b;
	Dm_c = U0 * H / Pe_c;
	Dm_d = U0 * H / Pe_d;

	A_a = tau - 0.5 - Dm_a / (cs*cs*dt);
	A_b = tau - 0.5 - Dm_b / (cs*cs*dt);
	A_c = tau - 0.5 - Dm_c / (cs*cs*dt);
	A_d = tau - 0.5 - Dm_d / (cs*cs*dt);

	printf("Nx=%d  Ny=%d\n", Nx ,Ny );
	printf("dx=%e tau=%e A_a=%e A_b=%e A_c=%e A_d=%e\n", dx, tau, A_a , A_b , A_c , A_d);

	k = 0;

	Init_Eq();
	while (k <=step)
	{
		k++;
		coll_BGK();
		Coll_MRT();
		Streaming();
		Boundary();
		if (k % step_out==0)
		{
            sprintf(matlab,"matlab%d.dat",k);
            sprintf(data,"data%d.dat",k);
            sprintf(data_a,"A%d.dat",k);
            sprintf(data_b,"B%d.dat",k);
            sprintf(data_c,"C%d.dat",k);
            sprintf(data_d,"D%d.dat",k);
			printf("k=%d  c=%e \n", k, c_a[Nx/4][Ny / 2]);
			data_output();
		}
	}
	data_output();

	printf("L_d=%d", L_d);

	system("pause");
}
